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ABSTRACT 

The statistical distribution of the radial pairwise peculiar velocity of galaxies 
is known to have an exponential form as implied by observations and explicitly 
shown in iV— body simulations. Here we calculate its statistical distribution 
function using the Zel'dovich approximation assuming that the primordial 
density fluctuations are Gaussian distributed. We show that the exponential 
distribution is realized as a transient phenomena on megaparsec scales in the 
standard cold-dark-matter model. 

Subject headings: cosmology: theory — galaxies: clustering — gravitation - 
large-scale structure of the Universe — methods: statistics and analytical 
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1. Introduction 

The pairwise velocity distribution function of galaxies plays an important 
role in translating observed clustering properties in redshift space to those in 
real space. Since Peebles (1976) proposed to adopt a distribution function of an 
exponential form, a number of observational analyses (Davis & Peebles 1983, 
Hale-Sutton et al. 1989, Fisher et al. 1994, Marzke et al. 1995) have confirmed it 
provides a better fit than other distributions. The results of iV-body simulations 
(Efstathiou et al. 1988, Fisher et al. 1994) have shown more directly that the 
distribution has indeed an exponential tail with a flatter peak. More recently, 
Zurek et al. (1994, hereafter ZQSW) have shown that the pairwise halo velocity 
distribution function has an skewed exponential shape with a sharp peak in 
their simulations of cold-dark-matter (CDM) models with a higher resolution. 

In spite of the simplicity of the functional shape, its origin had been a 
mystery for a long time. It was expected that the exponential shape arose as a 
result of highly nonlinear gravitational effects, because linear evolution preserves 
the initial shape of the velocity distribution which is usually assumed to be 
Gaussian. 

Recently Sheth (1996) and independently Diaferio and Geller (1996) 
proposed a model which explains the exponential shape of the pairwise velocity 
distribution function. They consider a universe composed by isothermal clumps 
whose mass function is given by Press-Schechter (1974) theory or observation. 
Assuming that velocity distribution of galaxies in each clump is Maxwellian 
whose dispersion is determined by its mass and that each pair of galaxies of 
interest is contained in a common clump, they have calculated the weighted 
sum of the pairwise velocity distribution function and shown that it has indeed 
an exponential tail. Diaferio and Geller (1996) have further discussed that if 
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they adopt a velocity distribution function of galaxies in a clump which is more 
highly peaked than Gaussian, they can reproduce the sharp exponential core of 
the final pairwise velocity distribution. 

Although their idea is interesting and attractive, its applicability is limited 
to a highly nonlinear regime with a small separation up to R <^ lh^ 1 Mpc 
(Diaferio & Geller 1996), because their approach is valid only when each pair 
belongs to the same nearly virialized clump whose typical size is around \hr x 
Mpc. On the other hand, the results of iV— body simulations clearly show that 
the exponential shape is kept for much larger separations up to ~ 5.5/i -1 Mpc 
(ZQSW) and for scales with a small value of the two-point correlation function 
£ ~ 0.1 (Efstathiou et al. 1988, Fisher et al. 1994). This motivates us to consider 
a different approach to explain the exponential shape on these scales; that is, 
not only a highly nonlinear effect but also a semi-nonlinear dynamics may play 
a role in realizing the exponential feature. 

In the present paper we calculate the pairwise velocity distribution function 
in the semi-nonlinear regime using a Lagrangian perturbation theory known as 
the Zel'dovich approximation (Zel'dovich 1970), assuming that the primordial 
density fluctuation is random Gaussian. In this approximation each mass 
element moves in a remarkably simple manner. In particular, peculiar velocity 
field remains proportional to the "shift" vector which is determined by the 
profile of initial density fluctuation. Hence, starting from Gaussian fluctuation, 
the statistical distribution of one-point peculiar velocity field remains Gaussian 
(Kofman et al. 1994). The pairwise velocity distribution of an arbitrary fixed 
separation, however, is obtained by an appropriate weighted sum of initial 
distribution function and can result in markedly non-Gaussian distributions in 
semi-nonlinear regime. As a result, an exponential shape is observed in the 
course of evolution as seen below. 
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The rest of the paper is organized as follows. In §2 we formulate the pairwise 
peculiar velocity distribution function in Zel'dovich approximation as a function 
of the initial relative velocity distribution assuming that it is Gaussian. Then 
in §3 we calculate it explicitly for the standard cold-dark-matter (CDM) model 
spectrum. We first trace its time evolution and then compare our results 
corresponding to the present epoch with N— body simulations done by ZQSW. 
Finally §4 is devoted to discussion and conclusion. 

2. Pairwise velocity distribution in Zel'dovich approximation 

In the Zel'dovich approximation (ZA), the Eulerian position, x(t), of a mass 
element is given by a function of its Lagrangian or initial coordinate q as, 

x(t) = q + D(t)p(q), (1) 

where p(q) is a shift vector and D(t) is a growth factor. The density contrast, 
5(x, t), reads 

-l 

-l~-D(t)V q -p(q), (2) 

where the latter approximate equality holds when S is small. D(t) is usually 
identified with the growth factor of the linear growing mode, and we do so 
here, too, but it can be an arbitrary linear combination of linear growing and 
decaying modes in principle. Then from equation ([]]) the peculiar coordinate 
velocity also has only the growing mode in the former case; 

v(t) = x(t) = D(t)p(q), (3) 

so it should be identified with the growing longitudinal velocity. 

Let us consider time evolution of elements A and B with the Lagrangian 
coordinates q A and q B = q A + r, respectively. The Eulerian position and 
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peculiar velocity at each point read 

x A (t) = x(q A , t) = q A + D(t)p{q A ) = q A + D(t)p A , 

x B (t) = x(q B , t) = q B + D(t)p(q B ) = q B + D(t)p B , (4) 

v A (t) = v(q A ,t) = D(t)p A , 

v B (t) = v(q B ,t) = D(t)p B . 

We investigate time evolution of their relative velocity, 

v AB (t) = v B (t) - v A (t) = D(t)(p B - p A ) 

= v\\(t) + v ± (t), (5) 

where v\\(t) and vj_(t) represent components parallel and perpendicular to 
r AB (t) = x B (t) - x A (t). They read 

f.y. v AB (t) ■ r AB {t) u ,r AB {t) 
Mt) = \r AB (tW rABit) = V ^VMt)> (6) 
v±(t) = V AB (t) - V\\(t), 

respectively, where r AB (t) and v\\(t) are given by 

r 2 AB (t) = + + (w V± <) ' (7) 

„ m 2 ( ^ -4- DV ABi \ (fl \ 

Here suffix i denotes quantities at the initial time t = ti when D{ti) = Di was 
negligibly small. 

Thus one can calculate the probability distribution function (PDF), 
P(V;R,t), of pairwise peculiar velocity v\\(t) = V with separation r AB (t) = R 
from the initial PDF, p(v v± xi , v± yi ] r), as 

P(V;R,t) = -^-^ J 4:7Tr 2 drdv\\ i dv ±X idv ± yip(v\\ i ,v ±xi ,v ± yi]r)5(R-r AB (t))5(V 

(9) 
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Here v± xi and v± yi are the two components of v±i vertical to each other. From 



the global isotropy the initial PDF depends on and v ±i = yjv\ jxi + v\ yi only. 
We find after some manipulation, 

2ttA 3 



Z7T 7J Z" 00 

P(y;R,t) = —+ rdrp(vl t ,vl t ;r), 

D 2 D Jr* 11 



(10) 



with 



*(r-°v 



and 



1 / £> x 2 



D 



(11) 
(12) 

(13) 



In order to evaluate the desired PDF let us specify the initial PDF, 
p(v*^, v* ±i ; r). Since we are dealing with only the longitudinal mode, in the 
initial or linear regime the peculiar velocity is related with S(x,t) or its Fourier 
transform, Sk(t), as 



3 ) 

2 



A J k 2 (2tt) 
(Peebles 1980). Assuming the initial density fluctuation is random Gaussian 
with the power spectrum Pi{k) defined by 
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we find the initial pairwise peculiar velocities are also Gaussian distributed with 
the two-point correlation functions (Gorski 1988), 
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One can also express Q\\(r) and Qj_{r) in terms of the integral of the initial 
two-point correlation function, £i(r), as 

Mr) 
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Q±(r) 
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(Peebles 1993, §21). We thus find 
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We can therefore obtain the desired PDF, P(V; R, t), through the integration in 
equation ([TI]) if we specify &(r) or Pi{k). The final expression reads 
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We readily see that it is independent of Z)«. Therefore we can deal with an 
extreme case peculiar velocity is set to be zero initially as in ZQSW without 
any difficulty Furthermore, since S»{r) and S±(r) are proportional to Df in the 
linear regime when the initial condition is set, equation (p2|) is also independent 

of A- 

Since Sn(r) and S±(r) approach the same value 

8tt 



Soo = yJ dkPi(k), 
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in the limit r — > oo, we can find an analytic expression of P(V; R,t) in the 
case r* is much larger than the typical correlation scale as 
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independent of R. Thus the PDF as calculated from ZA (eq. [p2] ) approaches 
Gaussian (eq. [p^] ) not only for a large separation, i£, but also in the cases V or 
D is large. The latter cases, however, correspond to a highly nonlinear regime 
where the Zel'dovich approximation is not reliable. 

Note that p{v*^v* Li \ r) in equation ([H]) has been normalized so that velocity 
integral of P(V;R,t) in equation (|22] ) results in the two-point correlation 
function of the evolved density field, £(.R, t), as 

/oo 
P(V;R,t)dV-l, (25) 
-oo 

(Bharadwaj 1996). 



3. Pairwise velocity distribution in CDM model 

3.1. CDM power spectrum 

In this section, we calculate the equation (|22"D with a specific model. We 
adopt the standard CDM model as a typical model of structure formation 
which has also been employed by ZQSW, namely, consider the Einstein de 
Sitter universe with the current Hubble parameter H$ = 50km/s/Mpc. We also 
normalize the scale factor a = 1 at present and set D = a. The initial power 
spectrum of density fluctuation is taken from Efstathiou et al. (1992) as 

m = — — ^ - T , (26) 

{l+ [ak + ((3k)l + ( 7 fc) 2 ] Y 

with a = 25.6 Mpc, (3 = 12 Mpc, 7 = 6.8 Mpc, and v = 1.13. B is the 
normalization factor which may be determined from the anisotropy of microwave 
background radiation as 
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where Q rms is the quadrupole fluctuation amplitude and T = 2.73K is the 
current temperature. The four-year COBE data gives Q rm s-ps — 15.3/xK 
(Gorski et al. 1996), while it ranges between Q rrns = 9.5/iK and 23.2/zK in the 
simulations of ZQSW. 

3.2. typical time evolution 

First we show how the distribution function of a fixed comoving separation, 
R, evolves in time. Figure 1(a) depicts evolution of P(V; R,t) for the comoving 
scale corresponding to R = 2 Mpc today at D — 0.05, 0.1, 0.2, 0.5, 1.0 and 
10.0. The horizontal velocity axis is normalized by the one point velocity 
dispersion proportional to DD and the numbers on the velocity axis apply 
for D = 1 only. Figure 1(b) represents time evolution for the comoving scale 
R = 6 Mpc today. These two figures imply that the negative tail develops much 
faster than the positive part from symmetric Gaussian distribution, properly 
reflecting the attractive nature of gravitational interaction. In the intermediate 
stage of the evolution around D = 0.5 — 1 for Fig. 1(a) and D ~ 1 for Fig. 
1(b), an exponential shape is prominent around the peak. We thus find that the 
exponential pairwise peculiar velocity distribution can be realized as a result of 
semi-nonlinear dynamics of the Lagrangian perturbation theory. 

In ZA, however, each mass element moves kinematically (eq. |l|) even after 
shell crossing, hence structures or correlations are destroyed as we extrapolate 
the evolution with it. As is seen in Figs. 1 the distribution becomes practically 
Gaussian again by the epoch D — 10 except around V = 0. In this highly 
nonlinear regime ZA is by no means valid. We thus see that non-Gaussian 
shape of the distribution function P(V;R,t) with a nearly exponential form 
near the peak is realized in this approximation as a transient phenomenon in 
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the semi-nonlinear regime. 

Finally we briefly comment on the evolution of the two-point correlation 



function of the density field £(R;t) obtained from equation ( j25j) . On the 
comoving R = 6 Mpc scale it evolves from to with the maximum value of 1.0 
at D = 1.2. ZA predicts smaller amplitude of £(R;t) than N— body simulations 
in the nonlinear regime (Schneider & Bartelmann 1995, WeiB et al. 1996). We 
return to this point below. 



3.3. Dependence on the separation 

Next we set D = 1 and analyze dependence on the separation R. We plot 
the new distribution function for the CDM spectrum normalized by COBE 
data, Q rm s = 15.3/zK, in Fig. 2(a) on scales R = 2, 6, and 10 Mpc. In Fig. 
2(b) the radial velocity distribution given by the linear theory is depicted 
for comparison. Even on the distance of 10 Mpc, these two distributions are 
significantly different from each other. 

As mentioned in the last subsection, with the original CDM power spectrum 
ZA predicts smaller magnitude of the correlation function in the nonlinear 
regime due to the fact that mass elements continue streaming even after shell 
crossing without getting bounded and virialized. To cure this problem it has 
been suggested in the literatures (Coles et al. 1993, see also Kofman et al. 1992) 
to adopt the "truncated" Zel'dovich approximation (TZA), in which small-scale 
power of the initial density fluctuation is phenomenologically suppressed using 
an appropriate window function such as a Gaussian filter and then apply ZA to 
the processed power spectrum, 

P t (k) = exp(-k 2 R})P i (k), (28) 

where Rf is the filtering length to be determined so that the results agree with 
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the predictions of iV— body simulations. Here we choose the filtering scale 
Rf = 2.2 Mpc so that the amplitude of the two-point correlation function at 
R = 6 Mpc becomes equal to that of the simulation (ZQSW). Figure 2(c) 
is the distribution function based on TZA where we find too small velocity 
dispersion. The filtering length we have adopted here is a modest one. In fact, 
for the purpose of optimizing the overall performance of ZA in reproducing 
the clustering properties of CDM iV— body simulations, an even larger value, 
Rf = 5.5 Mpc, has been suggested by Schneider and Bertelmann (1994). Then 
the magnitude of velocity dispersion would be even smaller. We thus find that 
TZA mimics properties of actual matter distribution by suppressing the velocity 
field and so it is not appropriate to analyze pairwise velocity field with this 
approximation. 

Next we compare our results with the numerical results of ZQSW who 
plotted frequency distributions of pairwise velocity for the pairs with separations 
1-2 Mpc, 5-6 Mpc, and 10-11 Mpc with the same condition as the figures 
of ZQSW with Qrms = 9.5/iK. In Fig. 3(a) we plot our result for R =1 Mpc 
(dotted line) and 2 Mpc (solid line) as well as the Gaussian tail distribution 
function Pa(V;t) defined in equation (|24]) by short-dashed line. Also depicted 
there is a long-dashed line which mimics the histogram of ZQSW near the peak 
where exponential (i.e. linear in the figures) fit is good. This line has been 
drawn from eye-fitting of Fig. 7(a) of ZQSW. Similarly in Figs. 3(b) and 3(c) 
results for R =5 and 6 Mpc and R = 10 and 11 Mpc are depicted together with 
approximate curves inferred from the N— body experiment of ZQSW. 

As is seen there we observe a systematic deviation to right in all the 
three scales. This is because gravitational infall is not completely taken into 
account in ZA. In fact the average infall velocities calculated by ZQSW are 
(Kium) = — 280km/s, — 355km/s, and — 276km/s for the pairs with R = 1 - 2 
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Mpc, 5-6 Mpc, and 10 - 11 Mpc, respectively, while ZA gives (Vza) — 14km/s, 
— 0.49km/s, and — 69km/s for R =2, 6, and 11 Mpc, respectively. If we correct 
this deviation of average infall by shifting the origin of the velocity axis with an 
appropriate amount, we find that the distribution function calculated with ZA 
reproduces the results of numerical simulation well especially near the peak in 
Figs. 3(b) and 3(c). 

This means that even if ZA fails to reproduce the average infall, it predicts 
pairwise velocity dispersion, ay(-R), fairly close to the corresponding numerical 
values. For an symmetric exponential distribution, the velocity dispersion 
appears in the PDF as 

V2\V\ 



P(V; R)dV oc exp 



dV. (29) 



a v (R) 

Since the PDF obtained from A— body simulations is asymmetric, it would 
be more appropriate to treat left- and right-hand-sides of the peak separately. 
Applying the fitting function of the type of equation (p5| ) independently on 
each side of the peak in Figs. 3(a)-3(c), we estimate the velocity dispersion 
on each side, oyi eftiright (-R), from the slope of each long-dashed line which is 
shown in the first column of Table 1 named Figure. Also shown there are the 
results based on ZA. The second column, ZA2, is obtained from the least-square 
fit between the theoretical curve (eq. fl2"2"f ) and equation (|29"D in the region 
P(V; R) < e~ 2 meixP(V] R) where the two sides are again treated independently. 
The third column, ZA1.5, is the same as ZA2 except for the narrower fitting 
region P(V; R) < e -1,5 maxP(V^; R). In both cases the central flatter region has 
been removed from fitting for R = 11 Mpc just as we have done so in drawing 
long-dashed curves in Fig. 3(c). 

As is seen there, ZA reproduces the asymmetric feature of the slope on both 
sides qualitatively well and on larger scale the agreement is even quantitative. 
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4. Discussion 

In the present paper we have formulated the pairwise peculiar velocity PDF 
with ZA starting from Gaussian fluctuations. Using the standard CDM power 
spectrum we have calculated the PDF on various scales at different epochs. 
As a result we have found that even if the one-point velocity distribution 
remains Gaussian, non-Gaussian feature develops prominently and we find 
an exponential shape on megaparsec scales before and around the epochs 
corresponding to the present. 

One may be tempted to calculate the relative velocity PDF of pairs of 
higher density regions or density peaks because the numerical data of ZQSW 
are those of halos where galaxies are likely formed. In ZA the density contrast 
is a function of dpi/dqj while velocity field is that of p iy and since they have no 
correlations in the isotropic Gaussian distribution, velocity PDF is independent 
of the configuration of the density field. 

Previous explanations on the exponential PDF (Sheth 1996, Diaferio & 
Geller 1996) are based on highly nonlinear dynamics of galactic systems. That is, 
the exponential feature appears as a result of superposition of Gaussian velocity 
distributions in clumps with various velocity dispersions. Hence although this 
model is suitable to account for a symmetric exponential distribution with 
vanishing mean net relative velocity, a more elaborate treatment is necessary in 
order to reproduce an asymmetric exponential distribution with its peak at a 
negative relative velocity as observed in numerical simulations. 

The PDF based on ZA naturally realizes an asymmetric exponential feature, 
although it also failed to predict the location of the peak correctly at D — 1. In 
fact, as discussed in the literatures (Coles et al. 1993, Schneider & Bertelman 
1994, WeiB et al. 1996) ZA does not successfully reproduce clustering properties 
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of density field in the regime D ~ 1, and it has been suggested to use the 
truncated spectrum to make the agreement better for D > 0.2 (Schneider 
& Bertelman 1994). On the other hand, as we have shown in Fig. 2(c), this 
TZA predicts too small amplitudes of the velocity dispersion. Thus neither the 
original ZA nor TZA reproduces both density distribution and pairwise velocity 
distribution simultaneously for the present universe. Hence it is not very certain 
to which extent the semi-nonlinear dynamics discussed here plays a role in 
realizing the actual exponential distribution in the present universe. In order 
to draw a more quantitative conclusion, it is desirable to compare our formula 
(eq.[[E|]) with the time evolution of N— body simulations of a CDM model, 
especially around the epoch D ~ 0.2 up to when the original ZA suffices for 
reproducing the power spectrum of density field. Unfortunately, we have been 
unable to do so due to the lack of published data with the sufficient informations 
about the initial condition, the epoch, and the separation in physical units. 

Our conclusion is therefore limited to a qualitative one. Nonetheless we 
believe it worth pointing out that simple kinematics as ZA results in an 
exponential shape of the pairwise peculiar velocity distribution. Indeed the 
pairwise PDF, which contains much more information on clustering in that its 
integral gives the two-point correlation function, is totally different from the 
one-point PDF. 
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FIGURE CAPTIONS 

Figs. 1 Evolution of the radial pairwise velocity distribution P(V;R,t) on 
the comoving separations R = 2 Mpc (a) and R = 6 Mpc (b) in terms 
of ZA with the CDM initial power spectrum with Q rms = 15.3/iK. Solid 
line represents the distribution at D = 0.05, dotted line at D = 0.1, 
short-dashed line at D = 0.2, long-dashed line at D = 0.5, dash-dotted 
line at D — 1 (corresponding to present epoch), and thick solid line at 
D = 10.0. To compare the statistical shape at different epochs, the velocity 
axis has been normalized by the one-point velocity dispersion proportional 
to DD and the numbers on the velocity axis apply for D = 1 only. 

Figs. 2 The distribution of radial pairwise distribution P(R; V, t) on the 

separations R = 2 Mpc (thin solid line), R = 6 Mpc (dotted line), R = 10 
Mpc (dashed line) and R = +oo (thick solid line) at present, D — 1. 
(a) represents the results based on ZA with the full CDM spectrum with 
Qrms = 15.3/iK, (b) is based on the linear theory with the same spectrum, 
(c) is obtained from ZA but with the truncated CDM power spectrum with 
the filtering length Rf = 2.2 Mpc. 

Figs. 3 Schematic comparison of our formula with the results of iV— body 
simulation of ZQSW for three typical separations at D — 1 with 
Qrms = 9.5/iK. In each figure, solid and dotted lines are the PDF obtained 
from ZA. The long-dashed line mimics the result of ZQSW in the regions 
an exponential fit is reasonable, (a) is for the pairs of the separation 
R = 1 - 2Mpc, (b) for R = 5 - 6Mpc and (c) for R = 10 - UMpc. The 
short-dashed line is the Gaussian distribution Pg(V; t) defined in equation 
(11)- 
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TABLE 1 

Velocity Dispersion Obtained from the Slope of the Exponential 

Fitting Function 

R (Mpc) Side Figure" (10 2 km/sec) ZA2 b (10 2 km/sec) ZA1.5 C (10 2 km/sec) 



2 left 8.3 5.7 5.9 

2 right 4.3 4.5 4.5 

6 left 8.5 7.0 8.1 

6 right 4.2 5.1 5.5 

11 left 7.0 7.0 8.2 

11 right 3.4 4.1 3.5 



a Based on the long-dashed lines in Figs. 3(a)-(c). b From the least-square 
fit between the theoretical curve (eq. f22|) and equation (|29"D in the region 
P(V; R) < e~ 2 maxP(V; R). c Sane as ZA2 but with a narrower fitting region, 
P(V;R) < e" L5 maxP(V A ; R). 
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